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Abstract 

The  concept  of  desingularization  in  three-dimensional  boundary  integral  computations  is 
reexamined.  The  boundary  integral  equation  is  desingularized  by  moving  the  singular  points 
away  from  the  boundary  and  outside  the  problem  domain.  This  allows  the  surface  integrals, 
which  become  nonsingnlar,  to  be  evaluated  by  simpler  techniques  and  speeds  the  computa¬ 
tion.  The  effects  of  the  distance  of  desingularization  on  the  solution  and  the  condition  of  the 
resulting  system  of  algebraic  equations  are  studied  for  both  the  direct  and  indirect  versions 
of  the  desingularized  boundary  integral  methods.  Computations  show  that  a  broad  range  of 
desingularization  distances  gives  accurate  solutions  with  significant  savings  in  the  computa¬ 
tion  time.  The  desingularization  distance  must  be  carefully  linked  to  the  mesh  size  to  avoid 
problems  with  uniqueness  and  ill-conditioning.  As  an  example,  the  desingularized  indirect 
approach  is  used  to  study  unsteady  nonlinear  three-dimensional  gravity  waves  generated  by 
a  moving  submerged  disturbance;  minimal  computational  difficulties  are  encountered  at  the 
truncated  boundary. 


1  Introduction 


Boundary  integral  equations  ( B IE )  provide  a  powerful  method  for  the  solution 
of  linear,  homogeneous  boundary  value  problems.  The  method  employs  a  funda¬ 
mental  solution,  which  satisfies  the  differential  equation  (and  possibly  part  of  the 
boundary  conditions),  to  reformulate  the  problem  as  an  integral  equation  on  the 
boundary.  In  conventional  BIE  formulations,  singularities  of  the  fundamental  solu¬ 
tion  are  placed  on  the  domain  boundarv.  This  requires  special  evaluation  of  sinnu-  - l- 
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lar  integrands,  which  can  result  in  costlv  numerical  calculations.  In  time-dependent  '  5 
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nonlinear  free  surface  problems,1-3  a  boundary  integral  problem  is  solved  at  each  J*‘J C. 
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time  step.  Since  most  of  the  computation  time  is  devoted  to  the  boundary  integral 
problem,  an  effective  solution  method  is  critical  in  the  time-marching  procedure. 

When  the  singularity  of  the  fundamental  solution  is  placed  away  from  the 
boundary  and  outside  the  domain  of  the  problem,  a  desingularized  boundary  inte¬ 
gral  equation  (DBIE)  is  obtained.  We  will  show  two  advantages  to  this  desingu- 
larization:  a  more  accurate  solution  can  be  obtained  for  a  given  truncation,  and  a 
numerical  quadrature  can  be  used  to  reduce  the  computational  time  to  obtain  the 
algebraic  system  representing  the  discretized  boundary  integral  problem.  There 
are  two  types  of  nonsingular  boundary  integral  formulations:  direct  and  indirect. 
In  the  direct  method,  Green’s  second  identity  is  used  to  derive  the  boundary  inte¬ 
gral  equation,  and  the  solution  of  the  problem  is  obtained  directly  by  solving  the 
boundary  integral  equation.  In  the  indirect  method,  a  boundary  integral  equa¬ 
tion  for  the  singularity  strength  is  formulated,  where  the  singularity  distribution 
is  outside  the  problem  domain.  (When  the  boundary  integral  formulation  is  sin¬ 
gular,  the  direct  and  indirect  method  can  be  shown  to  be  equivalent4.)  These  two 
methods  are  formulated  in  the  following  section. 

The  first  use  of  a  desingularized  method  is  the  classical  wTork  by  Von  Karman5 
which  determines  the  flow  about  axisymmetric  bodies  using  an  axial  source  dis¬ 
tribution.  The  strength  of  the  source  distribution  is  determined  by  the  kinematic 
boundary  condition  on  the  body  surface.  Kupradze6  proposes  locating  the  bound¬ 
ary  nodes  on  an  auxiliary  boundary  outside  the  problem  domain.  Heise7  studies 
some  numerical  properties  of  integral  equations  in  which  the  singular  points  are  on 
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an  auxiliary  boundary  outside  the  solution  domain  for  plane  elastostatic  problems. 
By  applying  Green’s  theorem  to  the  solution  and  a  simple  wave  source  with  the 
singular  point  lying  inside  the  body  (i.e.,  outside  the  problem  domain)  and  using  a 
bilinear  expansion  of  the  source,  Martin8  obtains  unique  solutions  of  the  so-called 
“null-field  equations  for  the  water  wave  radiation  problems.”  Han  and  Olson9  and 
Johnston  and  Fairweather10  use  an  adaptive  method  in  which  the  singularities  are 
located  outside  the  domain  and  allowed  to  move  as  part  of  the  solution  process. 
This  adaptive  method  requires  considerably  fewer  singularities  than  the  number  of 
boundary  nodes,  but  it  results  in  a  system  of  nonlinear  algebraic  equations  for  both 
the  strength  and  the  location  of  the  singularities.  All  these  studies  can  be  classi¬ 
fied  as  nonsingular  methods  and  show  a  considerable  reduction  in  the  computation 
time. 

The  most  complete  discussion  of  a  desingularized  technique  is  given  by  Webster11, 
who  uses  a  triangular  mesh  of  a  singularity  distribution  (a  simple  source,  in  this 
case)  placed  somewhat  inside  the  surface  of  an  arbitrary’,  three-dimensional  smooth 
body  to  find  the  external  potential.  The  integration  is  performed  analytically  for 
each  triangle  using  a  linear  distribution  of  singularity  strength.  These  integra¬ 
tions  require  evaluation  of  (logarithmic  and  arctangent)  transcendental  functions. 
From  the  numerical  results,  Webster  concludes  that  “Submergence  of  the  singular¬ 
ity  sheet  below  the  surface  of  the  body  appears  to  improve  greatly  the  accuracy, 
as  long  as  the  sheet  is  not  submerged  too  far."  This  indicates  that  for  a  certain 
discretization  of  the  body  surface,  one  may  obtain  a  more  accurate  solution  by 
the  nonsingular  formulation  than  by  the  singular  formulation.  Although  Webstar 
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did  not  attempt  it,  a  simple  numerical  quadrature  can  be  used  if  the  distance  of 
desingularization  is  sufficiently  large.  This  will  significantly  reduce  the  computa¬ 
tional  effort  required  to  obtain  the  algebraic  system  representing  the  discretized 
boundary  integral  problem. 

Although  nonsingular  formulations  of  the  indirect  method  have  grown  popular 
recently,  few  studies  on  the  desingularized  direct  method  have  been  reported,  es¬ 
pecially  for  three-dimensional  problems.  Schultz  and  Hong12  obtain  a  nonsingular 
direct  formulation  by  moving  the  singularity  in  the  Cauchy's  integral  away  from 
the  boundary  in  two-dimensional  potential  problems.  Their  results  also  show  the 
advantages  of  the  nonsingular  formulation.  They  also  use  an  overdetermined  sys¬ 
tem  combining  the  real  and  imaginary  parts  of  Cauchy’s  theorem  or  using  “extra" 
evaluation  points  away  from  the  boundary  contour.  This  overdetermined  system 
can  exhibit  higher-order  convergence  than  the  determined  system  from  the  real  or 
imaginary  part  of  Cauchy’s  theorem. 

Fewer  investigators  report  on  applications  with  boundaries  extending  to  infinity, 
especially  for  wave  problems.  Jensen,  Mi  and  Soding13  solve  the  steady  nonlinear 
ship  wave  problem  by  the  indirect  method  using  a  simple  source  distribution  above 
a  free  surface,  with  minimal  difficulties  at  the  truncated  boundary.  While  they  use 
an  upw’inding  technique  as  a  form  of  the  radiation  condition,  we  find  that  the 
unsteady  problem  can  be  described  by  a  technique  that  uses  no  special  treatment 
for  fixed  time  at  the  truncated  boundary,  as  long  as  a  desingularized  method  is 
used. 
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After  formulating  the  two  desingularized  methods  in  section  2.  we  then  examine 

the  effect  of  the  distance  of  desingularization  on  a  simple  problem  with  an  infinite 

/ 

plane  boundary.  The  convergence  of  the  method  with  respect  to  mesh  size  and  the 
computational  time  are  compared  to  find  “optimum”  desingularizatidn  distances. 
We  then  use  this  method  to  calculate  the  unsteady  nonlinear  waves  generated  by 
a  source-sink  pair  moving  below  a  free  surface. 

2  Desingularized  Boundary  Integral  Method 

The  Laplace  equation  is  the  governing  equation  in  a  domain  D  bounded  by  T.  We 
assume  that  a  Dirichlet  condition  is  imposed  on  part  of  the  boundary  Ti  and  a 
Neumann  condition  on  the  remaining  boundary 


> 

11 

o 

(in  D), 

(1) 

o 

II 

(on  Td), 

(2) 

d<f> 

!fo~x 

(on  Tn). 

(3) 

where  60  and  \  are  known  functions  and 

is  the  outward  normal 

to  the  boundary 

Tj.  If  the  boundary  extends  to  infinity,  a  far-field  condition  is  required. 

The  desingularized  boundary  integral  method  separates  the  integration  and 
control  (i.e.,  evaluation)  surfaces,  one  of  which  is  the  boundary  of  the  problem.  In 
the  direct  method,  the  boundary  of  the  problem  is  the  integration  surface,  while 
in  the  indirect  method,  the  boundary  is  the  control  surface. 
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2.1  Direct  method 


The  boundary  integral  equation  in  the  direct  method  is  derived  from  Green’s  sec¬ 
ond  identity: 

/  /  jD  (**  -  4D  -  / 1  (ogj  -  vg)  JT  =  0,  (4) 

which  holds  for  any  two  functions  with  second  derivatives  continuous  in  D  and  the 
boundary  I\  If  <j>  is  the  solution  and  rp  is  chosen  as 

i'  =  7  7  ~y~7  •  (5) 

\Xp 

with  its  singular  points  xp  outside  D  and  x,  on  I\  the  volume  integral  in  (4) 


becomes  zero.  We  then  have 


II ':(*%-*£)« -0-  <6) 

Applying  the  boundary  conditions  (2)  and  (3)  to  (6)  and  moving  the  known  quan¬ 
tities  to  the  right-hand  side  give: 

/  /  <*(?,)  j-  (^~ [~ )  dT  -  /  /  — --  d~~^dT 
J  Jr„  dn\\xv-xq\/  J  Jvd\xv-xq\  dn 

=  -/  /  fej-f.-  1  W-y/7  7—  ~  — -j \  <ir.  (7) 

^  vrd  dn  \|xp  -  x,|y  dyrn|xp-x?| 

The  kernels  are  nonsingular  since  xp  and  x,  never  coincide.  The  integral  equation 
(7)  is  solved  for  <j>  on  Tn  and  d4>/dn  on  f  d- 


2.2  Indirect  method 

The  indirect  method  forms  a  solution  by  integra'.ing  a  simple  source  distribution 
over  some  surface  fi  outside  the  problem  domain.  By  applying  the  boundary 
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conditions,  (2)  and  (3),  we  obtain  a  boundary  integral  equation  for  the  unknown 
strength  of  the  singularities,  cr(x,), 

\  I  <r(xa)---~-dn  =  4>0(xc)  (on  rd)  (8) 

J  \XC  X j| 

and 

//0^4(  j^by)  "-**•>  (onr»)'  (9> 

where  x,  is  the  integration  point  on  the  surface  fl  outside  D ,  xe  the  control  point 
on  T  and,  again,  ^  represents  the  derivative  normal  to  I\  Once  the  singularity 
strength  is  solved,  the  solution  for  <j>  can  be  determined.  The  term  l/|xc  —  x5|  can 
easily  be  replaced  by  other  higher-order  singularities  (dipoles,  etc.)  with  little  ad¬ 
ditional  computational  effort  since  the  integrals  are  nonsingular  and  are  evaluated 
numerically. 

2.3  Difficulties  introduced  by  desingularization 

Desingularization  results  in  a  Fredholm  integral  equation  of  the  first  kind  which  can 
lead  to  uniqueness  and  completeness  problems  of  the  solution  as  manifested  in  the 
ill-conditioning  of  the  resulting  algebraic  system  if  the  desingularization  distance 
is  not  properly  chosen.  Uniqueness  can  be  a  serious  problem  for  the  direct  method 
since  the  solution  of  the  algebraic  system  is  solution  of  the  problem.  However, 
we  can  accept  multiple  constructions  (different  values  of  a)  of  the  same  d>  with 
the  indirect  method.  Webster11  shows  that  the  nonsingular  formulations  are  not 
significantly  less  well  conditioned  than  singular  formulations  and  that  completeness 
is  assured  if  the  singularities  of  the  numerical  method  are  placed  closer  to  the 
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surface,  than  any  real  singularity.  He  also  discusses  strategies  for  choosing  the 
desingularization  distance  for  closed  bodies.  The  method  we  describe  here  uses  a 
desingularization  distance  that  is  related  to  the  local  mesh  size  (described  in  more 
detail  in  section  3).  As  the  meshes  become  finer,  the  singular  points  approach  the 
boundary.  In  the  limit,  the  nonsingular  formulation  is  consistent  with  the  singular 
formulation.  For  example,  the  desingularized  kernel  (although  never  singular)  can 
be  shown  to  have  pointwise  convergence  everywhere  except  at  the  singular  point. 
The  numerical  integration  error  will  converge  if  the  singularity  approaches  the 
boundary  at  a  sufficiently  slow  rate  as  the  mesh  is  refined  (as  shown  in  Section 
3).  Thus,  the  properties  of  the  singular  boundary  integral  equations  still  apply  for 
both  methods. 

Desingularization  increases  the  condition  number  of  the  resulting  linear  system. 
In  most  three-dimensional  computations,  the  large  number  of  unknowns  requires 
am  iterative  solution  of  the  linear  system.  As  the  condition  number  increases,  one 
can  expect  a  loss  in  accuracy  or  an  increase  in  the  number  of  iterations.  Then  there 
appears  to  be  an  “optimum"  desingularization  distance.  If  the  singularity  is  too  far 
from  the  boundary  the  linear  system  will  be  poorly  conditioned,  and  uniqueness 
and  completeness  problems  occur.  If  the  singularity  is  too  close  to  the  boundary, 
numerical  integration  of  the  singular  terms  is  suspect  and  the  solution  may  not  be 
as  accurate,  even  if  “exact"  integration  is  used.  Unfortunately,  little  guidance  in  the 
selection  of  the  desingularization  is  available  except  for  Webster’s11  discussion  for 
axisymmetric  bodies.  If  the  solution  is  sensitive  to  the  desingularization  distance, 
the  nonsingular  formulation  would  not  be  very  practical.  We  will  show  that  this 
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is  not  the  case  in  the  following  examples. 


3  Examples 

First,  we  test  the  numerical  performance  of  the  desingularized  boundan'  integral 
methods  on  a  problem  in  which  the  potential  is  generated  by  a  dipcle  below  an 
infinite  flat  plane  with  <f>  =  0.  This  simple  problem  has  an  exact  solution  formed  by 
the  dipole  and  its  image  about  the  flat  plane.  This  problem  represents  the  solution 
to  the  first  time  step  of  nonlinear  waves  caused  by  an  impulsive  disturbance  of 
a  dipole  under  water.  We  then  apply  the  method  to  time-dependent  nonlinear 
waves  caused  by  an  underwater  disturbance  (a  moving  source-sink  pair  in  our 
calculation). 

The  direct  and  indirect  methods  ar^  compared  for  the  simpler  example.  In  the 
direct  method,  the  boundary  is  divided  into  rectangles,  and  a  solution  is  sought 
at  the  nodes.  Since  the  kernels  are  nonsingular,  the  integrals  can  be  evaluated 
using  Gaussian  quadrature,  where  the  solution  is  assumed  to  be  bilinear  over  each 
rectangle.  In  the  indirect  method,  the  integrals  of  the  singularity  distribution 
(8)  and  (9)  may  be  replaced  by  a  summation  of  isolated  singularities  without  an 
apparent  loss  of  accuracy  when  the  desingularization  distance  is  appropriate.  This 
does  not  require  integration  and  mapping  in  numerical  calculations.  Therefore,  the 
computations  will  be  less  complex  and  time  consuming  than  in  the  direct  method. 

A  collocation  method  is  used  to  solve  both  formulations  of  the  boundary  integral 
equations:  In  the  direct  method  we  satisfy  the  Green  theorem  at  chosen  points. 
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while  in  the  indirect  method  we  satisfy  the  boundary  conditions  at  the  collocation 
points.  We  propose  that  the  distance  of  a  singular  point  from  the  boundary  be 
given  by 

Ld  =  ld(Dm)a ,  (10) 

where  ld  is  a  constant,  Dm  is  the  local  mesh  size  (usually  the  square  root  of  the 
local  mesh  area)  and  a  is  a  parameter  that  must  be  chosen  carefully. 

To  test  the  convergence  with  respect  to  mesh  ref  nement,  we  evaluate  the  poten¬ 
tial.  which  is  due  to  a  known  constant  normal  dipole  distribution  within  a  square, 
flat  surface  (  —  1  <  x  <  1,  —  1  <  y  <  1,  2  =  0),  at  a  point  with  a  distance  given 
by  (10)  above  the  center  of  the  surface.  The  surface  is  subdivided  into  a  square 
mesh  and  a  2  x  2  Gaussian  quadrature  is  used.  A  third-order  convergence  should 
be  found  assuming  the  integrand  is  nonsingular  and  is  independent  of  the  mesh. 
Although  moving  the  singular  points  makes  the  integrand  nonsingular,  Eq.  (10) 
makes  the  integrand  depend  on  the  mesh  size,  and  hence  third-order  convergence 
is  not  assured.  Fig.  1  shows  the  convergence  of  tha  numerical  integration  for  three 
values  of  a:  0,  1/2  and  1.  As  long  as  q  <  1.  third-order  convergence  is  recovered 
although  the  numerical  integration  error  increases  with  a.  On  the  other  hand,  a 
should  be  greater  them  0  for  the  uniqueness  and  completeness  properties  of  the 
solution  of  the  integral  equation.  Therefore,  an  appropriate  a  lies  in  between  0 
and  1.  In  subsequent  calculations,  we  choose  a  value  of  q  =  1/2. 

In  the  results  that  follow,  a  Generalized  Minimal  Residual  Method  (GMRES)14 
is  used  to  iteratively  solve  the  system  of  equations.  We  find  this  method  to  be 
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generally  twice  as  fast  as  a  standard  conjugate  gradient  routine  for  nonsymmetric 
matrices.  A  sufficiently  small  value  of  the  convergence  tolerance  is  chosen  so  that 
it  introduces  negligible  error  in  comparison  to  that  introduced  by  truncation.  The 
computations  are  carried  out  on  an  Apollo  DN10000  workstation  using  16  digit 
arithmetic. 

3.1  A  dipole  below  a  <5  =  0  infinite  flat  plane 

In  this  example,  a  Dirichlet  condition,  <£  =  0,  is  imposed  on  the  z  =  0  plane.  A 
dipole  of  unit  strength  is  located  at  x 0 ( 0 ,  0 ,  —  1 ) ,  and  the  direction  of  the  dipole 
coincides  with  the  i  axis.  The  normal  derivative  d<p/dn  is  sought  on  z  =  0. 

In  the  direct  method,  the  boundary  integral  equation  for  this  problem  is  derived 
from  Green’s  second  identity.  Since  the  surface  at  infinity  does  not  contribute,  we 
obtain 

I  f  ^-rdT  = - - — 3,  (11) 

J  hi  dn  \xp  -  x,|  \xp-x0\ 

where  xv  is  the  control  point  outside  ZT  and  xq  the  integration  point  on  The 
computational  domain  is  discretized  into  a  uniform  square  mesh  defined  by  A 
nodal  points  after  truncating  the  z  =  0  plane  (Tj)  at  x  =  ±FL> c  and  y  =  R ^  with 
a  symmetry  condition  on  y  —  0.  The  integrations  are  performed  using  a  bilinear 
distribution  of  d<f>/dn  and  using  a  2  x  2  Gaussian  quadrature.  The  control  points 
are  placed  directly  above  the  nodal  points  at  a  distance  Ld  =  /«,( Ax)1/J,  where  Ax  is 
the  mesh  spacing. 


In  the  indirect  method,  the  potential  is  formed  by 


N 


47T  \x  -  xJ  ~1  |x  -  X,  1 


(12) 


i=i 


where  x4,  are  the  source  points  above  2  =  0.  The  boundary  condition  <f>  =  0  on 
2  =  0  results  in 
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<7, 
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i=l  Fc  X«| 


4?r  (*?  +  Vc  +  (2c  +  l)2) 
where  xc  =  (xc,yc,zc)  is  the  control  point  on  z  =  0. 

The  exact  solution  is 

—  1  x 


3/2  ' 


(13) 


d-(x,y,r)  = 


+ 


41  (*’  +  +  (z  +  1)!)3/2  4,r  (*’  +  »>  +  (z  -  D2)3'” 

The  exact  solution  for  the  normal  velocity  on  the  2  =  0  plane  then  becomes 

_  _3 _ x 

dn  exact  2v  (x2  +  y2  +  1)5/2 


(14) 


(15) 


The  results  of  the  normal  velocity  are  compared  to  (15)  by  the  RMS  error,  defined 
bv 


=  Q  -(f; 

i>  |=1  \  071  car  act  Oil 

The  mesh  configuration  is  shown  in  Fig.  2. 


1 


N 


do , 


,<9o, 


(16) 


Fig.  3  shows  the  RMS  errors  for  this  example  using  the  direct  and  indirect 
methods  for  three  different  values  of  A’  while  varying  desingularization  distances. 
The  direct  method  using  a  2  x  2  Gaussian  quadrature  shows  a  rapid  drop  in  error 
near  zero  desingularization  followed  by  increasing  error  as  U  increases.  In  the 
indirect  method  (using  isolated  sources  x„  rigth  above  the  control  points  xc)  ,  the 
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solution  blows  up  for  U  =  0  because  the  integration  and  control  points  coincide. 
However,  as  Id  increases,  the  errors  decrease  rapidly  and  the  solutions  are  rather 
insensitive  to  the  variation  of  Id  for  quite  a  large  range.  Eventually,  Id  becomes  too 
large  for  the  given  truncation  to  represent  the  second  term  in  (14)  well,  and  the 
errors  start  to  increase. 

We  have  performed  some  limited  computations  where  the  control  and  nodal 
points  are  staggered.  These  staggered  computations  give  similar  results  for  opti¬ 
mal  desingularization,  but  give  significantly  better  results  than  the  nonstaggered 
algorithm  only  when  the  desingularization  is  very  small.  Since  staggering  strategies 
are  difficult  to  choose  and  only  marginal  improvement  is  obtained,  we  avoid  the 
staggered  method.  We  also  tried  using  a  surface  distribution  of  sources  above  the 
z  —  0  plane  in  the  indirect  method  for  this  problem.  This  is  similar  to  Webster’s 
method11  except  that  the  integrals  over  the  source  surface  are  performed  using  a 
2x2  Gaussian  quadrature.  As  seen  in  Fig.  4,  the  results  are  better  than  those 
of  isolated  sources  for  a  large  range  of  Id  (from  0.1  to  3.0).  This  can  be  expected 
since  the  use  of  isolated  sources  implies  using  a  rougher  quadrature  for  the  inte¬ 
grals  over  the  source  surface.  However,  we  found  that  the  condition  number  of 
the  resulting  system  of  equations  using  a  surface  distribution  of  sources  is  greater 
than  that  using  isolated  sources  by  an  order  of  100.  Since  it  does  not  improve 
the  optimum  accuracy  but  takes  considerably  longer  to  form  the  algebraic  system 
and  the  iterative  solution  will  take  longer  for  the  distributed  source,  we  find  that 
nonstaggered,  isolated  sources  are  prefered. 
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A  comparison  of  different  techniques  used  to  evaluate  the  influence  function 
for  the  direct  method  is  given  in  Fig.  5.  As  expected,  for  small  /<*  the  error  us¬ 
ing  the  Gaussian  quadrature  is  larger  than  the  error  using  Newman’s  analytic 
integration15  because  the  integration  error  dominates  the  accuracy  of  the  solution. 
As  U  increases,  the  Gaussian  quadrature  integrates  the  smoother  integrands  ac¬ 
curately,  and  the  results  merge  with  those  of  Newman’s  approach  for  all  lj.  It  is 
fortuitous  that  in  the  middle  range  of  U,  the  results  of  the  Gaussian  quadrature 
show  less  error  than  those  of  Newman’s  approach — the  discretization  and  numeri¬ 
cal  quadrature  errors  apparently  tend  to  cancel  each  other.  A  more  accurate  3x3 
Gaussian  quadrature  brings  the  results  closer  to  those  of  Newman’s  approach,  as 
expected.  Although  not  shown  here,  increasing  the  number  of  nodes  N  also  brings 
the  results  using  the  Gaussian  quadrature  closer  to  those  of  Newman’s  approach. 

When  using  the  direct  method,  the  singular  formulation  (U  =  0  in  Fig.  5)  gives 
the  least  error  when  analytic  integration  is  used.  Since  the  direct  and  indirect 
versions  are  equivalent  if  the  influence  function  is  evaluated  analytically  for  the 
singular  case,  we  expect  the  singular  indirect  method  to  be  equivalent  to  the  sin¬ 
gular  direct  method.  However,  desingularization  in  the  indirect  method  greatly 
reduces  error  as  can  be  seen  in  Fig.  3  and  Fig.  4.  Desingularization  is  apparently 
more  beneficial  in  the  indirect  method  because  the  solution  of  many  problems  can 
be  represented  accurately  by  a  finite  number  of  some  singularities  located  outside 
the  problem  domain.  In  this  example,  one  negative  image  dipole  above  the  z  =  0 
plane  results  in  the  exact  solution.  Desingularization  may  be  more  problem  (or 
geometry)  specific  for  the  indirect  method.  The  direct  method  seems  to  be  less 
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sensitive  to  desingularization. 


The  effect  of  desingularization  on  the  condition  number  of  the  system  of  linear 
equations  is  shown  in  Fig.  6  for  the  direct  and  indirect  methods.  The  condition 
number  increases  exponentially  with  ld.  However,  a  poorly  conditioned  system 
does  not  necessarily  imply  an  inaccurate  solution.  In  fact,  minimal  error  occurs 
around  ld  =  3.0  in  the  indirect  method,  where  the  condition  number  is  quite  large. 
The  increased  condition  number  is  likely  to  increase  the  number  of  iterations. 
Fig.  7  compares  the  ratios  of  CPU  time  by  GMRES  to  that  by  an  LU  algorithm 
for  the  solution  as  a  function  of  ld  for  N  =  231.  The  LU  algorithm  takes  about  3.5 
seconds  to  solve  the  system.  An  error  tolerance  e  for  the  least  squares  residual  of 
the  equations  needs  to  be  specified  when  GMRES  is  used.  We  find  that  e  =  10-5  is 
sufficiently  small  in  all  the  computations  performed.  With  the  use  of  this  tolerance, 
the  CPU  time  by  GMRES  is  less  than  that  by  the  LU  algorithm  for  <  1.5. 
Newman’s  approach  requires  about  22  seconds  to  set  up  the  matrix  while  the  2  x 
2  Gaussian  quadrature  requires  only  about  6.4  seconds,  for  a  savings  of  about  70 
percent.  However,  the  CPU  time  for  solving  the  system  varies  from  3.0  seconds  to 
4.0  seconds  as  ld  changes  from  0.8  to  3.0.  Even  in  the  worst  case.  ld  =  3.0,  we  gain 
a  total  reduction  in  CPU  time  of  about  60  percent  over  the  singular  formulation 
(with  ld  =  0  and  Newman’s  approach).  In  the  indirect  method,  the  CPU  time  for 
the  matrix  set  up  is  only  about  1.3  seconds,  and  the  CPU  for  the  solution  of  the 
system  varies  from  1.7  to  4.0  seconds.  The  total  CPU  time  is  reduced  by  about  SO 
percent  with  the  indirect  method  for  this  truncation.  When  the  desingularization 
distance  is  small  (e.g.,  ld  <  1.5),  a  larger  c  does  not  result  in  a  significant  difference 
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in  RMS  error,  see  Fig.  S.  However,  the  CPU  time  is  significantly  reduced,  see  Fig.  9. 


Fig.  10  shows  the  effect  of  truncating  the  infinite  plane.  The  computational 
domain  is  extended  by  adding  uniformly  sized  meshes.  Both  the  direct  and  indirect 
methods  converge  quadratically  with  respect  to  the  extent  of  the  computational 
domain  for  this  problem,  as  expected.  Both  methods  also  converge  algebraically 
(approximately  linearly)  with  respect  to  A',  as  can  be  seen  in  Fig.  11,  in  which  the 
computational  domain  remains  unchanged  while  the  mesh  size  varies.  The  mesh 
convergence  is  algebraic  for  all  (not  shown  here). 

Fig.  12  shows  the  error  distribution  along  the  x  axis  for  both  methods.  Larger 
oscillations  in  the  solution  by  the  direct  method  are  observed  at  the  edge  of  the 
computational  domain.  This  may  be  due  to  the  neglecting  of  the  contribution 
from  the  integrals  over  the  far-field  closure  in  the  direct  method,  which  very  likely 
results  in  a  strong  global  effect  on  the  accuracy  of  the  boundary  integral  equation 
itself.  In  the  indirect  method  no  integrals  are  required  over  the  boundary;  the 
boundary  condition  is  enforced  within  the  computational  domain,  and  the  singu¬ 
larities  individually  satisfy  the  far-field  condition.  Thus,  one  may  expect  the  effect 
due  to  the  failure  of  satisfying  the  boundary  condition  outside  the  computational 
domain  to  be  smaller  than  the  effect  due  to  the  neglecting  of  the  far-field  closure. 

3.2  Waves  generated  by  a  source-sink  pair  moving  below  a  free  surface 

For  an  irrotational,  incompressible  flow  in  am  ideal  fluid,  the  Laplace  equation 
is  the  governing  equation  for  the  velocity  potential  6.  On  the  free  surface,  the 
nonlinear  dynamic  and  kinematic  boundary  conditions  are 
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D4>  1 

^  =  -=/  +  5V^V* 


(17) 


and 

(18) 

where  xj  —  ( xj,yj,zj )  is  a  position  vector  of  the  fluid  particle  on  the  free  surface 
and  is  the  material  derivative  following  the  fluid  particle.  The  z  axis  is  taken  as 
positive  upwards.  Initially,  there  is  no  flow  and  the  free  surface  is  flat.  The  flow  is 
generated  by  the  motion  of  a  source-sink  pair  that  starts  from  rest.  The  speed  of 
the  disturbance  pair  is  quickly  brought  to  a  steady  value  by  a  smooth  function  of 
time  to  avoid  high-frequency  content.  The  problem  has  been  nondimensionalized 
by  taking  the  depth  of  the  disturbance  h  and  the  gravitational  acceleration  g  to 
be  unity. 

The  free  surface  boundary  conditions  are  integrated  with  respect  to  time  to 
update  the  position  of  the  fluid  particles  (the  nodes)  on  the  free  surface  and  their 
velocity  potential.  The  velocities  of  the  fluid  particles  on  the  free  surface  at  each 
time  step  are  determined  by  the  desingularized  boundary  integral  method.  The 
location  of  the  singularities  x,  are  moved  at  each  time  step  to  a  distance  Ld  away 
from  the  nodes  and  normal  to  the  surface  (the  normal  direction  is  approximately 
determined  assuming  that  it  is  perpendicular  to  the  two  vectors  defined  by  the 
diagonal  lines  of  the  four  adjacent  nodal  points).  Material  movement  of  the  nodes 
is  beneficial  for  this  problem  because  they  tend  to  naturally  cluster  near  crests 
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where  curvature  of  the  surface  and  velocity  gradients  are  large.  Since  nodes  cluster 
at  these  points,  the  desingularization  distance  also  decreases  in  a  beneficial  way  as 
prescribed  by  (10).  We  use  the  indirect  method  with  isolated  sources  because  of 
its  computational  advantages  and  simplicity  discussed  in  the  previous  sections.  An 
additional  advantage  of  the  indirect  method  for  this  problem  is  that  the  velocities 
can  be  obtained  directly  once  the  strength  of  the  singularity  distribution  is  known, 
while  the  direct  method  requires  numerical  differentiation  to  obtain  the  tangential 
velocities. 

The  potential  is  expressed  as  a  sum  of  1)  the  source-sink  disturbance  pair  at 
a  distance  h  below  the  undisturbed  free  surface,  2)  the  image  disturbance  above 
the  undisturbed  free  surface,  and  3)  a  sum  of  Ar  sources  of  unknown  strength  in 
an  array  at  distance  Ld  above  the  disturbed  free  surface.  The  far-field  condition 
is  better  satisfied  when  the  image  of  the  disturbance  is  used  in  the  construction  of 
the  solution.  The  “integral”  equation  for  the  unknowns  ax  is 


I3-/  Xtource 


.V 


O, 


-c(t) 


1=1  lJ7  x>' 

o  {i) 


\x]  -  \x}  -  x\ 


g(0 


(19) 


where  o{xs)  is  known  and  cr(i)  is  the  strength  of  the  disturbance.  After  the  cr,  are 
determined,  the  velocities  of  the  fluid  particles  on  the  free  surface  can  be  calculated. 
A  fourth-order  Runge-Kutta-Fehlberg  method  is  used  in  the  nonlinear  free  surface 
integration. 
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The  method  is  first  applied  to  waves  generated  by  a  sufficiently  small  distur¬ 
bance  such  that  linear  wave  theory  is  a  good  approximation.  We  compare  the 
results  of  the  present  method  using  fully  nonlinear  free  surface  conditions  to  an 
“exact”  solution  computed  from  a  time-dependent  Green  function  for  a  Kelvin 
wave  that  satisfies  the  linearized  free  surface  condition16.  Then  we  study  nonlin¬ 
ear  free  surface  waves  generated  by  a  stronger  disturbance.  In  both  cases,  the 
disturbance  (the  source  and  sink  pair  lying  along  a  horizontal  line)  moves  hor¬ 
izontally.  The  distance  between  the  source  and  sink  is  0.1.  The  pair  moves  at 
speed  V(t)  =  Fr(  1  —e~4t)  with  the  midpoint  between  the  source  and  sink  initially 
located  at  point  (5,0,  —  1).  Fr  is  the  Froude  number  defined  by  V0/y/gK ,  where  V0 
is  the  terminal  velocity  of  the  disturbance.  Ft  is  chosen  to  be  1  in  this  example. 
The  free  surface  initial  mesh  grid  extends  from  0  <  z  <  20  and  0  <  y  <  7.5  (with 
a  symmetry  plane  at  y  =  0)  and  is  initially  divided  into  an  40  x  15  mesh.  The 
nodes  in  the  x  direction  are  equally  spaced  while  the  distance  between  the  two 
adjacent  nodes  in  the  increasing  y  direction  increases  by  10  percent.  The  nodal 
points  serve  as  the  material  points.  The  time-marching  is  conducted  following  the 
nodal  points.  The  isolated  sources  are  placed  approximately  perpendicular  from 
the  nodal  points  at  a  distance  determined  by  (10),  where  the  local  mesh  size  Dm 
is  the  square  root  of  the  average  area  of  four  adjacent  meshes,  / d  is  1.0  and  a  is 
0.5.  The  magnitude  of  the  disturbance  for  both  the  source  and  the  sink  is  defined 
by  a(t)  =  a0(\  -  e_4f).  where  aQ  =  0.05  for  the  linear  case  and  oD  =  0.75  for  the 
the  nonlinear  case.  The  time  step  is  0.2  in  the  time-marching  procedure. 

Fig.  13  shows  the  wave  elevations  along  the  symmetry  plane  (y  =  0)  at  t  =  1.0 
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by  the  nonlinear  calculation  and  the  “exact”  linear  calculation  (see  King16)  with 
the  weaker  disturbance.  The  time  convolution  integrals  in  the  linear  calculation 
are  performed  numerically.  As  seen,  the  nonlinear  and  linear  results  agree  very 
well.  Independent  computations  using:  a)  a  smaller  computational  domain  (with 
the  same  mesh  spacing  w'ithin  0  <  i  <  15  and  0  <  y  <  7.5),  b)  finer  mesh 
grids  (80  x  15  and  40  x  30  with  the  same  computational  domain)  and  c)  doubling 
the  time  increment,  result  in  negligible  difference  for  the  nonlinear  calculation. 
This  indicates  that  even  for  the  small  disturbance  example  studied  here,  the  small 
differences  in  Fig.  13  are  primarily  due  to  nonlinear  effects. 

Fig.  14  shows  the  wave  elevations  along  the  symmetry  plane  at  t  =  1.0  for  the 
nonlinear  case.  The  differences  between  the  nonlinear  and  linear  results  in  this  case 
are  due  to  the  nonlinear  effect  of  the  free  surface  conditions.  The  figure  indicates 
that  nonlinear  effects  are  stronger  near  the  second  crest.  A  three-dimensional  view 
of  the  nonlinear  waves  at  t  =  1.8  is  given  in  Fig.  15. 

To  study  artificial  effects  (such  as  unexpected  numerical  wave  reflection)  caused 
by  far-field  truncation  on  longer  time  simulations,  the  weaker  disturbance  was  sim¬ 
ulated  on  the  two  differently  sized  computational  domains  mentioned  above  with 
identical  mesh  spacing.  Fig.  16  shows  two  marks  on  the  undisturbed  free  surface 
indicating  the  edges  of  the  two  computational  domains.  Very  small  differences  are 
seen  in  the  two  computations  before  the  wave  front  hits  the  edge  of  the  smaller 
domain.  Moreover,  even  after  the  wave  front  passes  this  edge,  significant  differ¬ 
ences  are  not  apparent  until  the  first  crest  passes  the  edge.  After  that,  the  error 
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starts  to  propagate  into  the  entire  domain.  This  differs  from  our  attempts  to  solve 
the  identical  problem  using  the  direct  method  (not  shown  here)  where  extrane¬ 
ous  reflections  were  noticed  almost  imediately  at  the  truncated  boundary  where  a 
radiation-type  boundary  condition  (in  this  case  4>  =  0)  had  to  be  imposed.  The 
small  differences  in  the  computations  of  Fig.  16  are  remarkable,  especially  consid¬ 
ering  that  the  indirect  method  using  material  nodal  points  does  not  require  any 
spatial  derivatives  and  hence  the  nodes  on  the  edge  of  the  computational  domain 
do  not  require  any  special  treatment  such  as  one-sided  derivatives,  not  to  mention 
the  radiation-type  boundary  condition.  It  appears  that  the  indirect  method  would 
allow  the  use  of  a  considerably  smaller  computational  domain. 

4  Conclusions 

The  nonsingular  formulation  significantly  reduces  the  time  required  to  compute 
the  matrix  of  influence  functions.  The  resulting  system  of  linear  equations  is  still 
adequately  well  conditioned  to  allow  efficient  iterative  solution  of  the  linear  alge¬ 
braic  system.  Accurate  solutions  can  be  obtained  by  the  desingularized  boundary 
method  for  a  large  range  of  desingularization  distances  on  the  order  of  the  mesh 
size  (near  /<*  =  1).  It  was  also  found  that  the  indirect  method  is  more  efficient  than 
the  direct  method.  It  is  easy  to  code  and  takes  less  computation  effort.  In  addi¬ 
tion.  the  indirect  method  appears  to  perform  better  in  problems  with  boundaries 
extending  to  infinity.  Nonlinear  wave  calculations  using  a  time-marching  proce¬ 
dure  were  greatly  facilitated  using  the  desingularized  boundary  integral  method 
with  remarkably  small  numerical  reflections  at  the  computational  boundaries. 
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4  Comparison  of  surface  distribution  and  isolated  sources  for  indirect  method 
(i?o o  =  6.667,  N  =  231,  ld  =  1.0) 

. Isolated  source; - Distributed  source 

5  Comparison  between  exact  and  Gaussian  quadrature 

(Rx  =  6.667.  .V  =  231.  ld  =  1.0); - exact 

- 2x2  Gaussian  quadrature; - 3x3  Gaussian  quadrature 

6  Effect  of  desingularization  on  the  condition  number  (Rx  =  6.667,  N  =  231) 

- Direct  method: . Indirect  method 

7  Effect  of  desingularization  on  computational  time  ( R^  =  6.667.  .V  =  231) 

— . . Direct  method; . Indirect  method 
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8  Effect  of  iterative  tolerance  on  error  (R x  =  6.667,  N  =  231,  Id  =  1  ;  Indirect  method) 

- e  =  10~2; . £  =  10"5 

9  Effect  of  iterative  tolerance  on  computational  time 
(R0 o  =  6.667,  N  =  231,  ld  =  1.0  ;  Indirect  method) 

- c  =  1G"2; - £  =  10"5 

10  Convergence  with  respect  to  truncation  of  boundary  at  infinity  (Ax  =  0.6667,  ld  =.  \) 

- Direct  method; . Indirect  method 

11  Convergence  with  respect  to  number  of  nodes  (Rx  =  6.667.  ld  =  1.0) 

- Direct  method; . Indirect  method 

12  Error  distribution  along  the  x  axis  (R0 0  =  6.667,  N  =  231,  Id  =  1.0,  y  =  0) 

- Direct  method: . Indirect  method 

13  Free  surface  elevation  along  the  symmetry  plane  at  t  =  1  (linear  case;  cr0  =  0.05) 

- Nonlinear  result; . Linear  result 

14  Free  surface  elevation  along  the  symmetry  plane  at  t  =  1  (nonlinear  case:  a0  —  0.75) 

- Nonlinear  result; . Linear  result 
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15  A  prospect  view  of  the  nonlinear  wave  at  t  =  1.8  (<70  =  0.75) 


16  Effect  of  computational  domain  size  on  wave  elevations  along  symmetry  plane 
{a0  =  0.05) 

- small  domain  (0  <  x  <  15); . large  domain  (0  <  x  <  20) 
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Direct  Method 

jp:  Control  Point 
xq:  Node  Point 


Indirect  Method 


xe:  Control  Point 
x Source  Point 
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Fig.  16 


